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Abstract 



1 Introduction 



We propose a new two-dimensional meshing algorithm called PINW able to 

generate meshes that accurately approximate the distance between any two domain 

points by paths composed only of cell edges. This technique is based on an extension 

of pinwheel tilings proposed by Radin and Conway. We prove that the algorithm 

produces triangles of bounded aspect ratio. This kind of mesh would be useful in 

| cohesive interface finite element modeling when the crack propagation path is an 

I> ' outcome of a simulation process. 

00 ' 

o 
t> 
o 

One of the most widely used techniques to simulate fracture is cohesive interface finite 
element modeling. In this kind of model, the area or volume under consideration is 
subdivided into bulk elements, which are typically triangles or quadrilaterals in 2D and 
tetrahedra or hexahedra in 3D. Next, interfacial elements, which are edge elements in 2D 
or surface elements in 3D, are placed between some or all pairs of adjacent bulk elements. 
The cohesive model prescribes a relationship relating traction to relative displacement on 
the interfacial elements. There is an abundance of literature that deals with the nature 
of this relationship, e.g., see and the references therein. A widely accepted modeling 
assumption is that the total energy to create the crack is proportional to its surface area 
(or length in 2D). In fact, the critical energy release rate G c per unit surface area or length 
of crack is often a parameter of the cohesive model. 

In a finite element model, the energy release rate is associated with surface area or 
length of interfacial elements composing the crack being modeled. If the discrepancy 
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between the "true" crack path (i.e., the path the crack would follow if it were not for 
the finite element constraint that the crack path must lie on predetermined interfacial 
elements) and the path of the simulated crack is large for certain paths, then nonphysical 
preferred crack directions can exist. In other words, the results of the simulation would 
depend upon how well the boundaries of the mesh cells are aligned along the true crack 
path. In this paper, we propose a meshing technique that approximates the true path 
with the path along mesh boundaries with high accuracy even though the true path is 
unknown to the mesh generation algorithm. In particular, the approximation has the 
property that for any crack path, the simulated and true crack path lengths converge to 
each other upon refining the mesh, which is a property not possessed by other simpler 
families of meshes. We call this algorithm the PINW mesh generator because it is based 
on an extension of the 1:2 pinwheel tiling described in the next section. 

In Section El we define "deviation ratio" and consider a simple experiment to test the 
properties of the 1:2 pinwheel mesh. The 1:2 pinwheel tiling seems to be too restricted to 
be useful for a general-purpose algorithm, so we explain how to generalize it to arbitrary 
triangles in Section |3J This generalization is the basis for our meshing algorithm PINW. 
In Section El we prove that our generalization still has the isoperimetric property. Then in 
Section El we describe the algorithm. The main new ingredient introduced in that section 
is a procedure to convert a tiling to a mesh. The aspect ratio of the resulting mesh is 
analyzed in Section 

The aspect ratio of the mesh is important for the cohesive fracture application because 
the bulk elements (e.g., triangles in 2D) are used to model a continuum mechanical theory 
such as linear elasticity. It is well-known (see, e.g., Theorem 4.4.4 of 0, in which aspect 
ratio is called "chunkiness" ) that poorly shaped elements can lead to substantial errors 
in the elasticity solution. 

2 Pinwheel tilings 

In this section, we provide a brief introduction to the properties of pinwheel tilings. Tilings 
are a covering of the euclidean 2-space £ 2 starting with a finite number of shapes called 
prototiles. The tilings are constructed by translated and rotated copies of the prototiles 
that intersect each other only along the boundaries. The tilings were proposed to model 
crystallographic structures in the physics community. 

The pinwheel tilings jH] are classified as aperiodic tilings. In S 2 this is equivalent 
to saying that no translation of the tiling leaves it invariant. The basic pinwheel tiling 
as developed by Radin and Conway has a hierarchical structure and is constructed by 
successive operations of subdivisions and expansions. 

Consider a right triangle with legs of length 1 and 2 referred to as the short and 
medium sides. The hypotenuse is thus of length a/5 and will be called the long edge. The 
vertices will be named similarly, that is, the small, medium and long vertices are opposite 
the corresponding sides. For brevity, we will call a right triangle with the ratio of its short 
to medium edge equal to 1/2 as a "1 : 2 right triangle" and the tiling formed by its copies 
as a "1 : 2 tiling." This single tile is subdivided into five triangles that are all congruent 
to each other as shown in Figure ^ 
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If one were to dilate the subdivision in Figure ^ by a factor of \/5 and then rotate 
and translate the resulting figure so that the dilated copy of C ended up coincident with 
the original tile P, then a larger subset of S 2 would now be tiled. The above subdivision 
scheme is then applied to each of the five copies of P, and then another dilation followed 
by rotation and translation is carried out. Continuing this process infinitely would lead to 
a tiling of the plane. Thus, in the case of the standard pinwheel tilings, P and Pr (where 
Pr denotes the reflection of P about the x-axis) form the set of fixed prototiles and the 
tiling uses translations and rotations of this set. 

For our purposes however, we will concentrate just on the subdivision step and omit 
the dilation, translation and rotation steps, leading to the "subdivision" pinwheel tiling 
in which the cell diameter tends to zero and the area of the plane covered by the mesh 
does not expand from step to step. This is because we are interested in generating a mesh 
with varying amounts of refinement for a fixed region rather than a mesh that ultimately 
covers S 2 . In the subdivision pinwheel tiling, one starts with a fixed 1 : 2 triangle and 
then repeatedly subdivides first the initial triangle and then each subtriangle into five 
congruent subtriangles using the above rule. 




-0.2 1 1 1 1 1 1 

0.5 1 1.5 2 2.5 

Figure 1: Basic pinwheel subdivision proposed by Radin. 

One can enumerate the rotation angles of the child triangles with respect to P and Pr 
as RqPr, R n +eP, ReP, RePn, R^+eP where Rg is rotation by 9 in the counterclockwise 
direction. For the standard 1 : 2 right triangle, 9 = arctan(l/2) and 9/tt in this case is 
irrational. The significance of this is as follows. As the number of subdivisions goes to 
infinity, so do the distinct orientations of the triangles. For example, suppose we keep 
track of the orientation of all triangles of type C with respect to the parent triangle in 
the subdivisions. As can be seen in Figure ^ the angle made by a triangle of type C 
with respect to the parent in the nth subdivision is nd. Since 9/tc is irrational, n9 will 
represent a different angle for each n. 

This presence of an infinite number of orientations leads to a special property known 
as the isoperimetric property For a tiling of S 2 , isoperimetry means that given an 
e > 0, there exists an R such that for any two points P and Q on the boundaries of the 
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triangles with | \P — Q\ \ > R, the shortest path from P to Q that uses only tile edges has 
length at most (1 + e)||P — Q\\. Here \\P — Q\\ denotes the Euclidean distance from P to 
Q, which will also be denoted as \PQ\. 

There is an analogous property for the subdivision pinwheel tiling. In this case, let 
P, Q be two points on the boundary of the initial triangle. Then for every e > 0, there 
exists an n such that after n recursive subdivisions of the initial triangle, the shortest 
path from P to Q using only triangle edges is at most | \P — Q||(l + e). This theorem can 
be generalized so that P and Q do not have to be on the boundary of the initial triangle 
but may be any two distinct points. 

The isoperimetric property is the reason that pinwheel tilings are attractive for cohe- 
sive interface modeling. Consider a finite region Q C S 2 tiled with an infinite sequence 
of pinwheel tilings M.\,M.2i ■ ■ ■ m which the triangles in M.i all have side lengths hi, 
2hi, \fEhi, and hi — > as i — > oo. Then for an arbitrary straight segment of length / 
connecting p G Q to q G £1, and for an arbitrary e > 0, there exists an I such that in each 
of the tilings JUj, Aii+i, . . ., there is a path from p to q using only mesh edges (except 
for initial and ending segments to connect p and q to the boundaries of the triangles that 
contain them) such that the length of the path is 1(1 + e). We will give a proof of this 
result in a more general setting in Section 

Since the above result holds for an arbitrary line segment, it also holds for any piecewise 
smooth curve or network of such curves. The reason is that a network of piecewise smooth 
curves can be approximated arbitrarily accurately with a network of line segments. Then 
each of the line segments can be approximated arbitrarily accurately with paths of the 
pinwheel tiling. 

Thus, when used for cohesive fracture, the pinwheel tiling has the property that all 
possible crack paths are approximated as accurately as desired (in terms of their length) 
by paths that use only mesh edges, as the mesh diameter tends to zero. As we shall see 
in the next section, more common mesh generation techniques do not have this property. 

3 A computational experiment 

In this section we carry out some simple experiments to quantify the isoperimetric prop- 
erty of the 1:2 tiling. Since our interest here is in meshes, we first explain how to convert 
the 1:2 pinwheel tiling to a mesh. It is apparent from Figure that the pinwheel tiling is 
almost a triangulation except for the hanging node bisecting the medium side of triangle 
E. We define a hanging node of a planar subdivision into triangles to be a point that is a 
vertex of one triangle but lies on the strict relative interior of an edge of another triangle. 

It is fairly simple to make the pinwheel tiling a mesh we divide every triangle 
into two by joining its medium vertex to the midpoint of its medium edge. In fact, it is 
not necessary to split all the triangles, and in our example we have obtained a mesh by 
splitting a certain subset of the tiles. This splitting is done only on the finest level of the 
pinwheel subdivision. 

Our computational experiment is as follows. Starting from a 1 : 2 rectangle, we divide 
it into two 1 : 2 triangles and then apply the pinwheel subdivision n times to each of the 
1 : 2 triangles. Thus, the final tiling has 2 • 5 n triangles. The resulting tiling of the original 
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Table 1: Direct computation of deviation ratios for the first five levels of pinwheel subdi- 
vision. 



rectangle is then converted to a mesh using the technique in the last paragraph. 

Given a tiling T of a domain Q, let Skel(T) be the 1-skeleton of T, that is, the union 
of all edges of all triangles, and let V(T) be the set of all vertices of T. Let I be a positive 
parameter chosen small enough so that Q contains a disk of diameter I. We propose to 
evaluate isoperimetric quality of the triangulation with the following quantity, which we 
refer to as the I -path deviation ratio: 



Here, distskei(T)( - > •) means shortest distance among paths restricted to Skel(T). The 
notation \\p — q\\n means the geodesic distance from p to q, i.e., the shortest path among 
paths lying in Q. Thus, this quantity measures the maximum ratio between the paths in 
the mesh versus geodesic paths. Clearly for any mesh T of any polygon, dev;(T) > 1. The 
pinwheel mesh of the 1:2 rectangle has the property that for any I G (0, 1), dev;(PT m ) — > 1 
as m — > oo where VT m is the pinwheel tiling of the 1 : 2 rectangle after m levels of 
refinement. 

Our experiment is to evaluate devi(VT m ) for V7[, . . . ,V%. The results are depicted 
in Table ^ The worst-case shortest path is shown in Figure 121 

In contrast, consider the meshes in Figure E3 The deviation ratios of these meshes have 
lower bounds greater than 1 irrespective of the number of subdivisions. In particular, the 
lower bound is y/2 ~ 1.414 for the mesh in Figure Ofa)- For the mesh that was used by Xu 
and Needleman ^U] (one of the first papers on cohesive finite element modeling), which 
is shown in Figure and is sometimes called a "cross-triangle quadrilateral" mesh, the 
worst case deviation ratio can be shown to be approximately equal to 1.082 in the limit 
as the mesh cell size tends to 0. 

4 Generalization of Pinwheel Tilings 

The 1 : 2 pinwheel tiling discussed up to now was extended to a tiling with an arbitrary 
right triangle and its reflection as a prototiles by Sadun [Hj. The small angle of the prototile 
determines the finiteness of the orientations and sizes of the tiles in the tilings that are 
discussed in [8 . We now describe our approach to extend the pinwheel subdivision to 
arbitrary (non-right) triangles. 



n devi(P7;) 



1 1.3416 

2 1.1948 

3 1.1843 

4 1.1264 

5 1.0831 



dev(T) 
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Figure 4: Our generalized pinwheel subdivision of an arbitary triangle into five subtrian- 
gles. 

First we propose a way of subdividing a general triangle and show that any number 
of subdivisions would produce triangles similar to a finite set of prototiles. Consider the 
triangle shown in Fig. HJ We denote the vertices by A, B and C in clockwise order and 
the included angles at these vertices by a, b and c respectively Assume also a < c. 
First, draw the segment CF such that F is a point on AB and /LFCB = a measured 
counterclockwise from CB. From F draw FD such that D is on AC and ZDFC = b 
measured clockwise from FC. From D draw E and G such that E is on AB and G is 
on CF and ZADE = b clockwise from DA and Z.GDC = a counterclockwise from DC. 
Thus, we have a subdivision of a general AABC into five triangles of which /, III and 
V are similar to the parent and the remaining two II and IV are similar to each other 
but not to the parent. Note that we required a < c to make this construction but we did 
not require any ordering on b. 

Theorem 1. The above procedure for subdivision produces triangles with angles belonging 
to the set A% = {a, b, c} or to the set A2 = {a, c — a, ir — c}. 

Proof. This is obvious by simply checking all the angles in Figure |U and using the fact 
that angles of a triangle sum to ir. □ 

Theorem 2. If the above subdivision procedure is used recursively on the subtriangles, 
then any triangle produced has angles either from A\ or A 2 . 

Proof. One checks that if we define a' = a, b' = c — a, d = tc — c then {a', d — a',71 — d} = 
{a,b,c}. □ 

For the rest of this paper, we say that a triangle with angles {a, b, c} (listed in this 
order) is conjugate to a triangle with angles {a, c — a, tc — c}. The point of Theorem |21 is 
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that conjugacy is a symmetric relationship. We remark that if the original triangle is a 
right triangle, i.e., c = tt/2, then this triangle is similar to its conjugate. This is the case 
considered by jH]. 

These two theorems imply a procedure for subdividing any initial triangle 7\ = AABC 
with angles a, b, c. Assume a < b < c. Apply the first subdivision rule to get five smaller 
triangles. Then, for the three similar to T 1; reapply the same rule recursively. For the two 
conjugates, apply the other rule. For the conjugate triangles, we do not necessarily have 
the order a < c — a < n — c, but we do not need that order. We need only the inequality 
a < ii — c, which must be true since a + b + c = n. 

This procedure runs into a difficulty when c ~ a (i.e., the initial triangle is close to 
equilateral) because in this case the conjugate triangle will have a bad aspect ratio. We 
get around this problem as follows. If c ~ a, then we first subdivide the initial triangle 
into three about its in-center, that is, we join the in-center to the vertices of the original 
triangle and form three subtriangles. We use a cutoff in our algorithm: if c — a is less than 
the cutoff, then the preliminary tripartition is carried out. The cutoff for c — a is chosen 
to optimize the smallest angle. In other words, a parent is divided about the in-center if 
the smallest angle prior to division is smaller than after the division. Here smallest angle 
happens to be the minimum of the angles in the two sets A\ and A 2 for a given set of 
angles {a, b, c} and can be shown to be rs 0.4 rad. Thus, we take the cutoff to be 0.4 rad. 

5 Isoperimetric property 

This section is devoted to showing the result that the generalization of the pinwheel tiling 
introduced in the previous section obeys an isoperimetric inequality. The analysis and 
proof technique in this section closely follow the proof from jjj. The following is the key 
lemma in the proof of isoperimetry. 

Lemma 1. Let triangle T = AABC be as above. Assume a/ir is an irrational number, 
where a is the angle of T at A. Let 9 G [0, 2tt) and e > be arbitrary. Then there is 
a refinement of T following the above rules that contains a triangle edge e such that the 
angle between e and the x-axis lies in the interval (9 — e, 9 + e). 

Furthermore, the length of e is at least ((a, b, c, e)L, where a, b, c are the angles ofT, e 
is as above, (() is a fixed positive-valued function, and L is the longest side-length ofT. 

Proof. Observe that Triangle III in the above subdivision is similar to the initial triangle 
T but is rotated by angle a. Call this triangle T'. If this triangle is subdivided by the same 
rule again, there will be another smaller copy of T, say T", rotated by 2a etc. The infinite 
sequence a, 2a, 3a, . . . taken mod 2tc is dense in the interval [0, 2tt) by the assumption 
that a/(27r) is irrational. Therefore, for some sufficiently fine mesh, there is an edge e of 
triangle in the interval (9 — e, 9 + e). 

For the second part of the lemma, observe that for any e > there is an n = n(e, a) 
such that every point in [0, 2ir] is distance (mod 2ix) at most e from at least one point in the 
set {a, 2a, 3a, . . . , na}. Therefore, one of T, T", . . . , described in the last paragraph 
will have the desired edge e. The longest side-length of T is L; the longest side-length 
of X" is q(a, 6, c)L, where q is some universal function (not depending on anything other 
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than a, b, c) derived from our construction. By similarity, the longest edge of has 
length q(a, b, c) 2 L. Thus, if we define ((a, b, c, e) = 7(0, b, c)q(a, b, c) n ^ e,a \ where 7(0, b, c) 
is the ratio of the shortest to longest side length of T, then the length of e is at least 
C(a, 6, c, e)L. This proves the second part of the lemma is also satisfied. □ 

For the first main theorem of this section, we need one more definition. We say that 
a generalized tiling T of a triangle T refines another generalized tiling 7" of T provided 
that for each tile r of T, either r appears in T or a subdivision of r appears in T . This 
definition implies that V(T) C V(T') and Skel(T) C Skel(T'). The first main theorem 
for this section is as follows. 

Theorem 3. Let T = AABC be a triangle with angles a, b, c such that a < c and a/n is 
irrational. Let T°,T X , ... be an infinite sequence of generalized tilings of T generated by 
the rules above. For each i, let yi be the maximum tile diameter in T\ We assume the 
sequence of tilings has the following two properties: (a) T t+1 refines T l , and (b) yi — > 
as i — > 00. Let P, Q be any two points on the boundary ofT. Then 

lim dist (P,Q) = \PQ\. 

i^ooSkcl(T 4 ) 

In other words, every straight-line path connecting two points (P, Q) on the boundary 
of T is approximated with arbitrary accuracy by a path of edges of the tiling. 

Proof. In order to prove this theorem, we require a simultaneous analysis of tilings of a 
conjugate triangle. Therefore, let us change notation so that the original triangle is 71, 
its conjugate is T 2 , and there are two sequences of tilings with the above two properties, 
namely, T^, , . . ., which are tilings of Ti, and 7^°, T^ 1 , . . ., which are tilings of 

Without loss of generality, let us further assume that 7^° for v e {1, 2} is simply {T u }, 
and that each subsequent TJ is obtained from T*~ l by splitting exactly one tile (so that 
7^ has exact 4i + 1 tiles). This assumption is without loss of generality because we can 
take our original given sequence T^ 1 ,^ 2 , etc., and insert all intermediate tilings (i.e., if 
the original was obtained from T^ 1 via 9 subdivision operations, then we can insert 
6 — 1 intermediate tilings in the sequence). If we prove that the limiting property holds 
for the augmented sequence, then it certainly also holds for the original sequence. 

We make the following preliminary observation about generalized pinwheel tilings. If 
S is any tiling of T v for v = 1 or 2 obtained from the above generalized subdivision rules, 
then there exists an % such that 7~* refines S. 

We use the following additional definitions and notation to prove the theorem. 

• Let dT v denote the boundary of T v , v — 1,2. Thus dTi and dT 2 are each unions of 
three segments. 

• Define X u = dT u x dT u , for v = 1, 2. 

• For any p = (P,Q) G X v , let f v (p,n) denote dist S kei(T™)(-P, Q)- 
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• For p = (P,Q) G X„, let g v (jp,n) = f„(p,n)/\\P - Q\\, v = 1,2. If we took the 
maximum of this quantity over choices of p, we would arrive at a quantity analogous 
to the "deviation ratio" introduced above in Section 01 Clearly g u (p,n) > 1 for all 
p, n. If P = Q, then define g u (p, n) to be 1. 

• Let F v {p) = inf n > f v {Vi n )- Note that f v {p,n) is a nonincreasing function of n 
(because every edge of is covered by edges of T™ +1 for all n), so this inf is also 
the limit of the sequence. 

• Let G u (p) = inf n > g v (p, n )- Clearly G u (p) > 1 for all p. For the same reason as 
above, G is also the limit of the sequence. 

The main theorem now reduces to showing that Gi(p) = 1 for all p G X\. Our proof 
technique requires us to claim more strongly that G u (p) = 1 for all p G X u and for 
v = 1 and 2. 

• Let s u be the length of the shortest edge of T u for v = 1,2. Let t u be the length 
of the shortest edge among the five triangles that result from one application of the 
splitting rule to T u , and let p u = s u /t u . 

• Let X' v C X v denote points p = (P, Q) such that P,Qe dT v and such that 
distal (P, Q) > tu- Note that X' v is a compact set under the norm specified above. 

The reason for introducing X' v is that G v is continuous on X' v as the following argument 
shows. Observe that for p = (P,Q) G X' u , G u (p) = F v {p)/\\P — Q\\. The function F v is 
continuous on all of X v because it is a metric. The denominator \\P— Q\\ is also continuous 
and bounded away from on X' vl hence G v is continuous on this set. (Once the theorem 
is proved, then it is established that G u is continuous on all of X U) but this is not so easy 
to prove at this stage of the argument.) 

The following lemma shows that it suffices to analyze X' v rather than all of X v . 

Lemma 2. For u — 1,2, 

sup{G u (p) : p G X u } < max(sup{Gi(p) : p G X[} } sup{G 2 (p) : p G X' 2 }). 

Proof. Choose an arbitrary p = (P, Q) G X v . This proof will show that 

G v (jp) < max(sup{GiQo) : p G X[}, sup{G 2 (p) : p G X' 2 }). 

Taking the supremum on the left will prove the result. If p G X' v then the result is 
immediate since that value appears in one of the two terms in the right-hand side. So 
assume for the rest of the proof that p G X v — X' u . If P, Q lie on the same side of T u , 
then the left-hand side is 1 because fv{p,n) = \\P — Q\\ for all n in this case since the 
boundaries of T v are covered by the edges of for each n. Since the right-hand side is 
greater than or equal to 1, the result follows immediately. 

The last case is that P, Q are on distinct sides (and in particular, are not vertices of 
T u ). In this case they must be less than distance t v of the same vertex by definition of 
X' v . For the rest of the proof of this lemma, consider only the v = 1 case since the v = 2 
case is similar. 
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First, suppose that P, Q are both within distance t\ of A, the vertex whose angle is 
a. Without loss of generality, P lies on AB and Q lies on AC. Consider the sequence of 
tiles Hq —Ti, Hi, H2, . . . such that Hi is the tile from TJ that contains vertex A. Each of 
these tiles is similar to T\. The diameter of the iT's tends to as i —>■ 00 by assumption. 
Thus, there is a K such that Hk contains both P and Q but Hk+i fails to contain one 
or both of P or Q. Let u be the length of the shortest side of Hk- We claim that either 
\A — P\\ > u/pi or \\A — Q\\ > uj p\. The reason is that if both ||A — P\\ < u/p\ and 
1 1 A — Q || < u/pi then P, Q would both lie on the boundaries of H K+ i since the side lengths 
of Hk+i are all at least u/pi by definition of p\. This would contradict the choice of K. 

As mentioned above, Hk is similar to Ti, and the constant of proportionality is u/s\. 
Note that Hk could be either a dilation of Ti with no reflection or a dilation of Ti with 
a reflection. Assume the former case since the latter is similar. There exist P and Q 
lying on sides AB, AC of T\ whose positions with respect to AB, AC are proportional to 
the positions of P, Q with respect to the two sides of H K - Since either \\A — P\\ > uj p\ 
or \A — Q\ > u/px and the scaling factor between Hk and Ti is u/s\, this means that 
at least one of P, Q is distance from A greater than or equal to (u / p\) / (u / S\) = ti. 
Hence p = (P,Q) G X[. For an arbitrary n > 0, consider the tiling of Hk that is 
obtained by shrinking T x n by a factor of u/s\ and translating it so that it lies on top of 
Hk- The shortest path between P and Q in this tiling is fi((P ',Q),n)u/ ' s\ by scaling. 
Also, there is an n' such that the portion of Ti™ lying in Hk is strictly a refinement of TJ 1 
by the observation made at the beginning of the proof. The distance between P and Q 
in this tiling is /i((P, Q), n'), and since refines TJ 1 , f x {{P, Q),n') < f 1 ((P,Q),n)u/s 1 . 
Note that u/s\ = \\P — Q\\/\\P — Q\\ by similarity, so the previous inequality implies 
gi((P, Q),n') < gi((P, Q),n). Take the infimum over all n of both sides to conclude that 
Giip) < Gi(p), thus establishing the lemma in this case. 

In case that P, Q are both within distance t% from vertex B whose angle is b, the 
lemma follows by the same argument since the triangles containing B in all subdivisions 
of Ti are similar to T\. 

The last case is that P, Q are both within distance t\ of vertex C whose angle is c. 
Say, e.g., that P lies on AC and Q on BC. In this case, the argument is slightly more 
complicated since there are two triangles containing C in the next level of subdivision. Let 
CF be the segment that is the common boundary to the two triangles of the next level of 
subdivision that meet vertex C. (Refer to Fig. UJ) Let R be the point where segment PQ 
crosses edge CF. Then the argument above shows that the infimum over n of the distance 
between P and R using edges from Ti n is less than or equal to \PR\ ■ G2{P,R), where 
(P, R) G X' 2 . Similarly, the infimum over n of the distance between R and Q using edges 
from Ti n is less than or equal to \RQ\-G\(R, Q), where (R, Q) G X[. Therefore, T\(P, Q) < 
\\PQ\\ ■ max(sup{Gi(p) : p G X[}, sup{G2(j>) : p G X' 2 }) so the result follows. □ 

Finally, we conclude the proof of the main theorem by showing that sup{G v (p) : p G 
X' u } — 1 for v — 1, 2. To this end, choose v (either 1 or 2) so that sup{Gr„(p) : p G X' u } > 
sup{Cr3_j,(p) : p G X' 3 _ u }. Without loss of generality, say u — 1 is chosen. 

Since X[ is a compact set and G\ is continuous on this set, there exists a p = (P, Q) in 
X[ that maximizes G\(p). If P, Q lie on the same side of Ti, then Gi(p) = 1 so the proof 
is finished. Else let the corresponding pair of points where the supremum is achieved be 
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Q 



Figure 5: Construction for the proof of Theorem El 



p* = (P*,Q*) and H* be the line segment joining them. Assume (for a contradiction) 
that G^p*) = S > 1. 

Choose N\ large enough so that there exists a tile T' in T^ 1 such that T" n H* has 
positive length and is contained in the middle third H*. Let the longest edge of T' be u\. 
For reasons to be explained below, we also choose N\ large enough so that 

Ui < (VS-l)C(a,b, c, aiccos(l/VS)) 

where £() is the function defined by Lemma H Continue splitting until we reach split 
number N 2 > N\ so that within T^ 2 , there exists a tile T" in T^^ 2 lying inside in T' that 
has an edge making an angle 9 with H*, with cos^ > 1/y/S. This is possible by Lemma ^ 
See Figure El Let L2 be the length of this edge. By the second part of the lemma, we 
may assume L 2 > C( a ; b, c, r/)iti, where 77 = arccos(l/\/5')- 

Let Vi, V2 be the endpoints of this edge with Vi being the vertex near P* and V2 near 
Q*. Observe that P*Vi, ViV 2 and V 2 Q* make up a three-segment path from P* to Q*. 
The length of this path is L x + L 2 + L 3 , where L x = \P*V 1 \ and L 3 = |V 2 Q*|. (We have 
already defined L 2 = | V1V2 1 -) Let /1, l 2 , and l 3 be the lengths of the projections of P*Vi, 
\iV 2 , V 2 Q* respectively onto H*. Because V\ lies within T' while H* crosses through T', 
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the distance from V\ to H* is at most U\. hence 



Ll < jl\+u\=hyjl+ulll{ 

< h(l+ul/ll) 



- ll+ h 

hence 

«? 3w? . . 

(The factor of 3 arises because \\P* — V\\\ > \\P* — Q*\\/3 as assumed earlier.) Similarly, 

L 3 - l 3 < (3) 
3 _ |P*Q*| v ; 

Next, because cos# > -4= where 8 is the angle between H* and V\ V2, we have V / S'/2 > -^2- 
Therefore, 

,SZ 2 - L 2 > (v 7 ^ - l)L 2 . (4) 

Next, note that Fi((P*, Q*)) < 5Lx + L 2 + SX3 thanks to the existence of the three- 
edge path P*V{V<iQ* ■ The reasoning is as follows. From P* to V\ there is a straight-line 
path of length L\. This path cuts through a finite list of triangles, say <fi triangles, within 
the tiling T^ 2 , since P* and V\ both lie on triangle edges of this tiling. Let the individual 
segments within these triangles be of length pi, . . . ,p^>. By construction, these quantities 
sum to L\. Then by further refinement, we can find paths within the tiling with lengths 
less than or arbitrarily close to piS, P2S, etc. since S is the factor that is the maximum 
amount longer that an edge path in refinements of either 7\ or T 2 can be versus the 
straight-line path. So the infimum of the lengths of these paths added up is at most SL±. 
The same reasoning accounts for the term SL 3 . Finally, the edge V\ Vi is length L 2 and 
is already in the tiling. 

Use ©, ©, © and the equality \P*Q*\ = h + h + h to bound Sin + I 2 + SI 3 : 

F l ((P*,Q*)) < SL^SIz + Ii 

= S{U - h) + S(I 3 - h) - (Sl 2 - L 2 ) + S(h + l 2 + l 3 ) 

< &ul-^—-(^/S-l)I 2 + S\P*Q*\. (5) 
Multiply (Q) by U\ on both sides and use the fact that I 2 > ((a, b, c, r))ui to obtain 

u\ < iVS ' ~ S 1)L2 \P"Q"\- (6) 

Substituting ©in © shows that Fi((P*,Q*)) < S\P*Q*\. But this is a contradiction, 
because the hypothesis of this analysis was that G\((P*,Q*)) = S, i.e., Fi((P*,Q*)) = 
S\P*Q*\. □ 
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The preceding theorem has the drawback that it pertains only to paths starting and 
ending on the boundary of the root triangle. For isoperimetry, we would like to generalize 
the result to paths with arbitrary interior P and Q. Since the nodes of the pinwheel tiling 
are dense in the interior (in the limit as the mesh size is refined), the following theorem 
provides a suitable generalization and will be taken as our definition of the isoperimetric 
property. 

Theorem 4. Let T°, T 1 , ... be a sequence of generalized pinwheel tilings ofT (satisfying 
a < c and a/ir is irrational as in the previous theorem) such that the maximum cell 
diameter tends to zero and such that T l+1 refines T l for all i = 0, 1, 2, . . .. Let P, Q be 
any pair of distinct points lying on Skel(T n ) for some n. Then 

lim dist (P,Q) = IIP -Q II. 

m-»oo Skel(T m ) 
m>n 



Proof. Consider the segment PQ lying in T. Let e > be given. Make a list U±, . . . , U r 
of tiles in T n traversed by this segment. Since PQ crosses Ui, define P%Qi to be Ui D PQ. 
Observe that Pi,Qi both lie on the boundary of Ui. By the preceding theorem, after a 
sufficient number of further subdivisions (say s), there exists a path in Skel(T n+s ) between 
Pi and Qi of length |Pj<5i|(l + e). This choice of s depends on i, so take the maximum 
such value of s (maximum over alH = 1, . . . , r). Then there is a path in Skel(T n+s ) from 
P to Q of length at most 

|PiQi|(l + e) + |P 2 Q 2 |(1 + e) + • • • + \P r Q r \(l + e), 
i.e., at most \PQ\(1 + e). □ 



6 Meshing an arbitrary region 

In this section we present our algorithm PINW to mesh a region Q with arbitrary polygonal 
boundary. A summary of PINW appears in Figure El The steps in this summary are 
described in more detail in the remainder of this section. The current version of PINW 
is 1.0 and has been coded in Matlab. An example output from this algorithm is shown in 
Fig.0 

We first start with a coarse triangulation of the domain. We use the Triangle package 
developed by J. Shewchuk, which uses Delaunay triangulation. The triangles produced 
have bounded aspect ratio. The second preliminary step, as mentioned in Section HJ 
locates triangles too close to equilateral and splits them at their in-center. 

A third preliminary step is to identify and split triangles whose smallest angle a 
is a rational multiple of n. (As noted above, the proof of isoperimetry requires that 
a/n be irrational.) In principle, this test could be conducted exactly using number- 
theoretic methods since the coordinates of the vertices of each triangle, being floating- 
points numbers, are rational numbers and can be treated with integer algorithms by 
clearing common denominators. Modifying a triangle in which a is a rational multiple 
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Algorithm PINW 1.0 

1. Generate a mesh for Q with bounded aspect ratio using Triangle. 

2. Split triangles too close to equilateral at their in-centers. 

3. Split triangles whose smallest angle is a rational multiple of n at a point 
near the in-center. 

4. Let the set of triangles obtained after steps 1-3 be called %. 

5. Initialize a heap containing triangles that need splitting. The triangles are 
ordered so that the one whose minimum altitude is maximum is at the top 
of the heap. Initially the heap contains all triangles from %. 

6. Repeatedly remove a triangle from the heap and split it into five, until the 
size of the top element of the heap is sufficiently small according to the 
user's specification. 

7. Let % be the set of tiles including those in % and all their descendants 
obtained by subdivision. Let Tf G % be the set of leaf tiles. 

8. Loop over all tiles in % starting from the coarsest to determine the value 
of big(e) for each edge e of any tile. 

9. For each big edge (i.e., each edge in the image of the "big" operator), select 
one side as moving and the other as staying. Sort the list of nodes lying on 
the staying side of each such edge. 

10. Loop over tiles in % starting from the coarsest excluding Tf. For each 
such tile T and for each of its vertices D,E,F as labeled in Figured let 
e be the maximal big edge containing the particular vertex. If this vertex 
D, E or F is on the moving side of e and is very close to a vertex v' on 
the staying side, then displace it to coincide with v' and apply the induced 
affine transformation to subtriangles of T. 

11. Apply Delaunay triangulation to each distorted, subdivided leaf tile. (The 
distortion of the leaf tiles is due to the affine transformations in the previous 
step. The subdivision of the edges is due to the presence of hanging nodes.) 
The collection of triangles output from this step is a simplicial mesh of Q. 



Figure 6: Overview of the steps of the PINW algorithm. 
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Figure 7: The coarse mesh for this figure had three triangles. The final mesh after 
pinwheel tiling, collapse-node operations and Delaunay triangulation is shown. 

of 7i is trivial in principle because any small random perturbation of a node of such a 
triangle will lead to an angle that is not a rational multiple of n with probability 1. 

In practice, this exact test and solution are both undesirable. For practical use of the 
algorithm, we would like to avoid the case when a is close to a rational multiple of n of the 
form mix/n where n is a small integer. The reason is that the presence of such a triangle 
in which a/n is irrational but is close to m/n implies that, although the isoperimetric 
property is asymptotically valid, the available angles will be badly distributed (clustered 
around multiples of n/n) for modest levels of refinement. 

Therefore, a more practical heuristic is to check each smallest angle against a finite list 
of the form rrnr/n, where m, n range over a pre-selected set of small integers. If a triangle's 
smallest angle comes too close to a member of this list, then the triangle is either split 
into three using a point near its in-center or is perturbed. (The exact in-center obviously 
should not be used since this would replace each angle by half its previous value, and 
hence still close to a rational multiple of tt.) This step has not been implemented in 
the current version of our code PINW 1.0 because we are still seeking the best practical 
heuristic. (Indeed, in Figure one coarse triangle is close to a 45-degree right triangle, 
and hence one part of the subdivision exhibits a shortage of possible directions.) 

Let Tq be the list of triangles that are produced by these preliminary steps. Thus, the 
triangles in Tq form a simplicial triangulation of the input set Q. We call these triangles 
the root tiles. The generalized pinwheel subdivision is then performed on the the root 
tiles to obtain a refined tiling. The procedure to refine the mesh used in PINW 1.0 is 



16 



based on a simple heap |Xj. The heap is initialized with all triangles in %, which are 
ordered in the heap according to length of the minimum altitude. The main loop for the 
subdivision is to remove the top member of the heap (i.e., the unsubdivided tile with the 
largest value of minimum altitude) and replace it with its five children. The procedure 
terminates when the top triangle in the heap is smaller than the user-specified mesh size 
requirement. 

Note that during the subdivision procedure, the angles a, b, c in Figure |U are assigned 
to smallest, middle and largest angles respectively for tiles similar to root tiles. For the 
conjugate tiles, angles a, b, c are assigned according to the conjugacy relationship. In other 
words, if the angles of the root tile are a', b', d in that order, then the angles a, b, c in the 
conjugate tile are assigned in the order a = a', b = c' — a', and c = n — d . This ensures 
that the conjugate of the conjugate is again similar to the root tile. 

From this description, it is apparent that PINW 1.0 supports a single global user- 
specified mesh size requirement. For many applications of mesh generation, it is useful to 
have a finer mesh in one part of the domain versus another. This can also be implemented 
in the framework of generalized pinwheel subdivision but is not available in PINW 1.0. In 
addition, several aspects of our analysis that follows below would have to be generalized 
to cover graded meshes. 

Once the subdivision procedure is complete, the resulting tiling must be converted to 
a simplicial mesh. For the 1:2 pinwheel triangulation, this step is quite straightforward 
as mentioned in Sectional In the generalized case, however, it is much more complicated 
and involves several steps that we shall now describe. 

Let %, be the list of all tiles in the hierarchy: it includes the tiles in T and all 
their descendants from the subdivision procedure. The tiles in % naturally have a forest 
structure associated with them in which the forest roots are root tiles. Let leaf tile denote 
a triangle in % that is not further subdivided during the generalized pinwheel subdivision 
phase. Let 7} be the set of leaf tiles. 

The first step in converting the tiling to a mesh is to identify for each edge e of each 
tile T G % the edge that we denote big(e). This is defined to be the edge d of a triangle 
T' higher up in the subdivision hierarchy (i.e., T is derived from T 1 via a sequence of zero 
or more subdivision operations) such that e C d, and such that d is maximal with this 
property (i.e., there is no other ancestor of T with an edge e" that strictly contains d). 

For each triangle in T , big(e) = e. For some other tile T with an edge e, it is a 
straightforward matter based on a checking a finite number of cases whether big(e) = e 
or big(e) ^ e. In the latter case, big(e) can be determined from the immediate parent of 
T (assuming big(e) is already tabulated for the the parent's edges). Thus, it is possible 
to determine big(e) for each edge of each tile in % with a constant number of operations 
per tile. 

Next, for each "big" edge e (that is, an edge such that big(e) = e), identify a moving 
and staying side. This choice can be quite arbitrary, except for two stipulations. An 
edge e adjacent on the exterior boundary of fl should have its inside labeled staying (i.e., 
no tiles lie on its moving side). An edge in correspondence with CF in Figure 0] (every 
big edge generated during the subdivision procedure is in correspondence with either DE, 
DF, CF or DG) should have the side facing vertex B labeled as moving. We now identify 
all the nodes on the staying side of e and sort them in order of occurrence on the edge. 
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(a) 



(b) 



Figure 8: Example of a collapse-node operation is shown. A node on one side of a "big" 
edge that lies within the tolerance of a node on the other side is moved and merged with 
the nearby node on the other side. 

This sorted list is saved for the next phase of the algorithm. 

In the next phase, we loop over triangles in % — Tf starting from the coarsest and 
perform collapse-node operations on each. Let T be a tile in % — Tf. Let the four vertices 
of T introduced when it is subdivided be labeled D, E, F, G as in Figure |U We perform 
no operation for G since it is on the staying side of edge CF. The maximal big edge 
containing D is big(AC); call this b(D). The maximal big edge containing E and F is 
big(AB); call this b(E) and also call it b(F). Let v be one of D, E, F. We check whether 
v is on the moving side of b(v ). If it is on the staying side, then no further operation is 
performed. If it is on the moving side, then we find the vertex v' taken from the staying 
side of b(v) that is closest to v. This v' can be found efficiently using binary search on 
the precomputed sorted lists. If ||t> — v'\\ < 5, we collapse nodes v and v' by displacing v 
to v'. Here 6 is a tolerance discussed more below. 

This displacement induces uniquely determined affine transformations on triangles 
contained in T as follows. If v is the vertex labeled D in Figure then there is a 
unique affine transformation on AADF that leaves A and F fixed and moves D to v'. 
A second affine transformation of ACDF leaves F and C fixed and moves D to v'. If 
v is the vertex labeled E, then there are unique affine transformations determined for 
AADE and ADEF. Finally, if v is the vertex labeled F, then there are transformations 
for each of ADEF, ACDF and ABCF. The algorithm applies all the relevant affine 
transformations caused by motion of the node. Note that the affine transformations agree 
on the boundaries between these triangles, so there is no consistency issue regarding which 
transformation to apply. These transformations move the triangle, including every node 
at deeper levels of the hierarchy contained in it. This concludes the description of the 
collapse-node operation. See Figure |H] for an illustration of this operation. 

Note that a single tolerance S is used to determine motion. The theoretical value for 5 
is given by (jHJ) below. We will verify later that this value of 8 is sufficiently small so that 
two important properties hold: 
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Property 1 of 5: If a vertex v' is the target of a collapse-node move, then it should be 
uniquely determined, i.e., there should not be two vertices v' and v" on the staying 
side of b(v) that are both within distance 5 of v. 

Property 2 of 5: No two vertices on the moving side of big(e) for any e should be 
collapsed to the same node on the staying side. 

In a future extension of PINW to handle graded meshes, presumably the value of 5 would 
not be a single global value. 

The affine transformations described above have the property that all of the segments 
illustrated in Figure |U remain straight (collinear) segments after the transformation. It 
is apparent that each collapse-node operation could cause many nodes to move. We will 
say that the one node v that is displaced to match v' is directly displaced. The nodes 
moved by virtue of an affine transformation induced by moving v are said to be indirectly 
displaced. 

A collapse-node operation, once executed, cannot be undone by future collapse-node 
operations. The reason is that v is never moved again. It is never moved again directly 
since it can be moved directly only when the tile T = AABC that created it is processed. 
It can also never be moved again indirectly since there is no tile in lower levels of the 
hierarchy that contains it except as a corner vertex, and corner vertices of a triangle T' 
are not moved when T" is processed. Similarly, v' can never be moved again. The reason 
is that v' is never moved directly (since it is on the staying side of its big edge). Any 
transformation that might move v' indirectly takes place at a level of the hierarchy higher 
than the level of T. 

We carry out all available collapse-node operations for all triangles in the order de- 
scribed. Once all collapse-node operations are complete, we are left with the collection of 
distorted leaf tiles, each of which may have one or more hanging nodes. These hanging 
nodes are collinear with the endpoints of the edges on which they hang because, as noted 
above, we do not disturb any collinearity relationships with collapse-node operations. The 
hanging nodes are all at least 8 apart from the corners and from each other. 

For each of these distorted tiles, we compute its Delaunay triangulation (including the 
hanging nodes). The collection of all of these Delaunay triangles forms a simplicial mesh 
that is the final output of PINW. 

The running time of PINW is analyzed as follows. Let n be the number of leaf tiles. 
Then the total number of tiles is 0(n), as is the total number of vertices and edges. 
The heap insertions and deletions require 0{n\ogn) total operations. Sorting all the lists 
associated with big edges requires 0{n\ogn) operations. Looking up a vertex in a sorted 
list requires O(logn) operations for binary search, hence all of the lookups to see if a node 
should be collapsed require O(nlogn) operations. 

The recursive application of affine transformations requires 0(nd) operations since 
each vertex is transformed at most d times, where d is the maximum depth of the forest 
associated with T*. We claim d = O(logn). It follows from Lemmas El and El in the next 
section that the minimum altitude of a triangle at depth k lies between a§C k and a\D k , 
where cxq, cti are lower and upper bounds on the minimum altitudes among root tiles, D 
is an absolute constant and C is a scalar depending on the worst aspect ratio among root 
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tiles. This means that a leaf tile can be at most a factor log Dj log C (asymptotically) 
deeper in the forest than any other leaf tile. Thus, all leaves have depth O(logn). 

Finally, the Delaunay triangulation operations in the last step of the algorithm also 
require O(nlogn) operations total. Overall, we see that PINW requires O(nlogn) oper- 
ations. 



7 Analysis of aspect ratio 

In this section we analyze the aspect ratios of triangles produced by PINW, showing that 
they are bounded above by a number that depends only on the sharpest angle in the 
original polygon Q. Before this analysis, we first explain how to select the parameter 
S described in the last section. The parameter 5 depends on the minimum altitude of 
leaf tiles as will be apparent from the theory developed here. Let minalt(T) denote the 
minimum altitude of triangle T. 

Lemma 3. Let T be a triangle and let a = minalt(T). Then T can be enclosed between 
two parallel lines at distance a apart. Conversely, ifT can be enclosed between two parallel 
lines at distance a apart, then minalt(T) < a. 

Proof. The first part of the lemma is quite trivial: draw a line through the longest side 
length of T and a parallel line through the opposite vertex. These lines are distance a 
apart. The argument for the converse is as follows. Without loss of generality, let the two 
lines be parallel to the x-axis. Let the vertices of T be numbered v±,v 2 , v 3 such that v 1 is 
closest to the bottom line (i.e., has minimal y-coordinate among the three vertices) and 
t>2 to the top line. By reflecting if necessary, assume also that the ^-coordinate of v\ is 
less than or equal to the rc-coordinate of v 2 . Now draw the line v ±v 2 , which is a transverse 
to the two parallel lines. If v 3 is below (to the right) this line, then it is easy to see that 
the entire triangle T may be rotated clockwise about v\ until Viv 3 becomes horizontal, 
and during this whole rotation, all three vertices remain between the lines. On the other 
hand, if v 3 is above (to the left) of the line, then rotate T clockwise about v 2 until v 2 v 3 
becomes horizontal. 

Once T has been reoriented so that one of its edges is horizontal, the claim is trivial 
since the altitude to the horizontal edge is a vertical line segment and hence must have 
length no more than a. □ 

Corollary 1. lfT\,T 2 are two triangles such thatT\ C T 2 , then minalt(Ti) < minalt(T 2 ). 



Proof. Draw the two parallel lines for T 2 as in the previous theorem; clearly T x also lies 
between them. □ 

Lemma 4. Let T be a triangle with vertices v±,v 2 ,v 3 . Let T' be the triangle with vertices 
v 'i, v 2,V3- Let A be the unique affine transformation that carries T to T' . Let I be an 
arbitrary line segment. Then 

lengthM(Z)) 

l - d ' a < - l L g t H \) - 1+d ' a < 7 > 
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where d = dist(vi,v[) and a is the altitude of v\ with respect to V2V3. 

Proof. Without loss of generality, assume a is 1 and d is replaced by p = d/a. Furthermore, 
without loss of generality, let T be positioned so that its v 2 v 3 edge is a subsegment of the 
x-axis and Vi lies on the y- axis (hence v\ = (0, 1) by the previous assumption). With these 
assumptions, the afhne tranformation A in the question becomes a linear transformation 
(i.e., no additive term) since the x-axis (and the origin in particular) is invariant. The 
transformation maps maps (1, 0) T to (1, 0) T and (0, 1) T = v\ to v[. Let r = v[ — V\, and 
let r be written as (a,/3) T so that a 2 + f3 2 = p 2 . Then A corresponds to the matrix 



The minimum and maximum distortion of a line segment under a linear transformation is 
governed by the minimum and maximum singular values of the transformation. Thus, the 
question now hinges on the two singular values of A. Notice that A may be regarded as a 
perturbation of the identity matrix, which has two singular values equal to 1. Therefore, 
by Corollary 8.6.2 of the largest singular value of A is at most 1 +p, and the smallest 
singular value is at least 1 — p. These values are attainable by taking a = and (3 = 
±p. □ 

Lemma 5. Consider the generalized pinwheel subdivision illustrated in Figure ||o/ a 
triangle T. Assume that b > min(.4, a) and c — a > min(.4,a). Then letting T' be any 
one of the five subtriangles, we have minalt(T') < 0.9725 minalt(T). 

Remark 1. The assumptions are valid for all tiles produced by PINW. For root tiles, we 
have ordered the angles a < b < c, and we know c — a > .4 because of the preliminary 
step of splitting near-equilateral triangles. For conjugates of root triangles, say a = a 1 , 
b = d — a', and c = tt — d where a', b', d are the angles of a root tile, we know also 
c — a = 7r — d — a' = b' > a and that b = d — a' > .4. 

Remark 2. The factor 0.9725 is due to our proof technique and appears to be an 
overestimate. A search over a fairly dense grid of possible angles satisfying the hypotheses 
of the theorem indicates that the true bound is closer to 0.918. 

Proof. We start by observing that |i?C|/|AB| = sin a/ sine by the law of sines. We know 
that either c > 2a or c > a + .4. Furthermore, we know that either ir — c > 2a or 
it — c > a + .4 since 7r — c = a + b. Now consider two cases. Case 1 is that c > tt/2. Define 
c = it — c, so that sin c = sin c, c > min(2a, a + .4) and c < tt/2. On the interval [0, tt/2], 
the sine is increasing. In the subcase that c > 2a, we have sine > sin 2a = 2 sin a cos a > 
V^sina. The last inequality follows because a < it /A (by assumptions that c < tt/2 and 
c > 2a) so cos a > l/y/2. The other subcase is that c > a + .4. Since sine is concave and 
increasing on [0, tt/2] the worst case (maximum value) for sin a/ sin c is when a = tt/2 — .4 
and c = tt/2, so sin a/ sin c < .922. Thus, \BC\/\AB\ < .922. The other case is c < tt/2. 
This case is handled by the same argument, except using c in place of c. 



Observe that subtriangle V, which is denoted Ty, is similar to T except scaled by a 
factor \BC\/\AB\. This proves minalt(T v ) < .922 minalt(T). 
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Next, by similarity, \BF\/\BC\ = \BC\/\AB\ hence \BF\/\AB\ < .849. This means 
\AF\/\AB\ > .151. Since AADF is isosceles, \AD\/\AF\ > .5 and \DF\/\AF\ > .5 hence 
\AD\/\AB\ > .075 and \DF\/\AB\ > .075. Next, by the law of sines applied to ACDF, 
\DC\/\DF\ = sin6/ sin(c — a). We now take three cases: either b < .4, b G [.4,7r/2], or 
b > it/2. In the first case b > a since the assumption in the lemma is b > min(.4, a). 
Also, c — a > it/ 2 since c — a + 2a + b = tt and 2a + b < 1.2 by the assumption for 
this case, so c — a > tt — 1.2. This means sin(c — a) = sin(7r — (c — a)) = sin(2a + b), 
with 2a + b < tt/2. Next, 2a + b < 3b so sin(c — a) = sin(2a + 6) < sin(3o) (since 
2a + b < 3b and sine is increasing on [0, 7r/2]), so sinfe/ sin(c — a) > sinfe/ sin(36). Now 
sin(36) = sin 6(3 cos 2 b — sin 2 6). Thus, sin 6/ sin(c — a) > 1/(3 cos 2 b — sin 2 6). Since 
sin 2 is increasing while cos 2 is decreasing on [0, .4], the minimum value of this fraction 
is when b = 0, so sin6/ sin(c — a) > 1/3. In the second case, b G [.4,7r/2]. This means 
sin b G [sin(.4), 1] and in particular, sin b > .389 so sin 6/ sin(c — a) > .389. The last case 
is b > tt/2, which implies sin(o) = sin(7r — 6) = sin(c + a). So the quantity to analyze is 
sin(c + a)/ sin(c — a). Since the angles in the numerator and denominator are both less 
than tt/2 (because b > tt/2) and c + a > c — a, we conclude that sino/ sin(c — a) > 1. 

Thus, in all cases, we conclude that \DC\/\DF\ > 1/3. This means that \DC\/\AB\ > 
.025. 

Next, observe that \EF\/\DF\ = sin(c— a)/ sine by the law of sines applied to ADEF. 
Again, we take three cases. If c — a < .4 (and hence c — a > a i.e., c > 2a, i.e., a < c/2), 
then 2a — a < .4 i.e, a < .4 so c < a + .4 < .8. Since all these angles are in [0, 7r/2], 
sin(c — a) > sin(c — c/2) = sin(c/2) and sin(c/2)/ sin(c) > 1/2 in this range by the 
convexity of sin. The second case is c — a G [.4, 7r/2]. In this case, sin(c — a) > .389 so 
sin(c — a)/ sine > .389. The last case is c — a > tt/2. This implies that sin(c — a) = 
sin(7r — c + a) = sin(2a + 6). The denominator becomes sin(c) = sin(7r — c) = sin(a + 6). 
So we are analyzing sin(2a + 6)/ sin(a + 6), which exceeds 1 since all the angles in question 
are in [0,tt/2]. Thus, in all cases, \EF\/\DF\ > .389 so \EF\/\AB\ > .0292. 

Now we have enough inequalities to analyze minalt(Ti) where T\ denotes AADE. 
Observe that this triangle is similar to T. Its longest side is either AE or AD (but not 
DE, since a < c). If its longest side is AD, then we see that \ AD\ = \AC\-\DC\ < \AC\- 
.0251^^1 < \AB\ - .025|AB| = .9725|AB|. (Here we used the fact that \AB\ > \AC\, 
which follows from the hypothesis of this case that \AD\ > \AE\ plus the similarity of 
AADE to AABC) Thus, AADE is similar to AABC but is a factor .9725 or less scaled 
down. 

The other case is when the longest side of ADE is AE. In this case, \AE\ = \AB\ — 
\EB\ < \AB\ - \EF\ < (1 - .0292)|AB| < (1 - .0292)|AC|. The inequality \AB\ < \AC\ 
follows from the assumption that \AD\ < \AE\ and similarity. Thus, ADE is similar 
to ABC but is scaled down by factor less than .9708. This concludes the analysis of 
minalt(Ti). This same analysis applies to Tin, since Tm is congruent to T\ (because 
AADF is isosceles). 

Next, observe that |AE|/|AD| = sino/sinc by the law of sines applied to AADE. 
Again, we take three cases. If b < A (and hence b > a), then c > tt/2 as above so 
sine = sin(7r — c) = sin (a + 6), so the quantity to analyze is sin 6/ sin (a + 6). Using 
analysis like before, including steps like sin(a + 6) < sin(26) = 2 sin b cos b < 2 sin 6, we 
conclude that sin 6/ sine > 1/2 in this case. If b G [.4, n/2], then we conclude again that 
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sin bj sine > .389. Finally, if b > tt/2, then sin b = sin(7r — b) = sin(c + a), so again 
sin b/ sine > 1. Thus, in all cases, \AE\/\AD\ > .389. Since \AD\/\AB\ > .075, we 
conclude that \AE\/\AB\ > .029. 

Next, we analyze Tn, that is, ADEF. Observe that this triangle is similar ACFA, 
and minalt(AC-Fy4) < minalt(Av4I?C) by Corollary d The corresponding side to AF 
is FE. We have \FE\ = \AF\ - \AE\ < \AF\ - .029|AB| < \AF\ - .029\AF\. Thus, 
minalt(Tn) < 0.971 minalt(T). 

Finally, we analyze Tjy, which is ACGD. This triangle is also similar to ACFA. We 
showed earlier that \CD\ < .963|AC|, and |CD|/|AC| is the ratio of similarity between 
these triangles. □ 

The following lemma is like the previous one except with an inequality in the opposite 
direction. 

Lemma 6. Consider the generalized pinwheel subdivision illustrated in Figure of a 
triangle T. Assume that b > min(.4, a) and c — a > min(.4, a). Then letting T' be any 
one of the five subtriangles, we have minalt(T') > pminalt(T), where for subtriangles I, 
II, III, IV, p > 0.0044 and for subtriangle V, p > sin a. 

Remark. The factor 0.0044 is due to our proof technique and appears to be an under- 
estimate. A search over a fairly dense grid of possible angles satisfying the hypotheses of 
the theorem indicates that the true bound is closer to 0.125. 

Proof. Again, we consider the five subtriangles and reuse some of the inequalities in the 
preceding proof. Starting with Ti, which is similar to AABC, recall that \ AD\ > .075|AB|, 
hence by similarity, mina\t( A AD E) > .075 mmalt(AABC). The same bound applies to 
Tni, which is congruent to Tj. 

For Ty, by the law of sines |i?C|/|AB| = sin a/ sine > sin a. Thus, by similarity, 
minalt(Tv) > minalt(T) • sin a. 

Next, consider triangle AACF. In the previous proof we showed that |AF|/|AB| > 
.151, which means that if AACF were dilated by 1/.151, it would completely cover 
AABC. Therefore, by Corollary d minalt(AACF) > .151 minalt(AABC). Next, we 
showed that \EF\/\AB\ > .0292 so \EF\/\AF\ > .0292. Since T n is similar to AACF, 
we conclude that minalt(Tn) > .0292 minalt(AACF) > .0292- .151 minalt(AABC). Note 
that .0292 • 0.151 > 0.0044. 

Finally, to analyze Tiv, we need to develop new inequalities. Recall we have already 
shown that \DF\/\AB\ > .075. Since ADGF is similar to AACB, this implies \DG\ > 
.075|AC|. Meanwhile, we know \AC\ > 0.5|y4F| since \AC\ is not the shortest side of 
AACF (because it is opposite an angle of size a + b, which is greater than the angle at 
A of size a). Thus, \DG\ > .037] AF|. This means by similarity of Tjy to AACF that 
minalt(T IV ) > .037minalt(AACF) > .037 • .151 minalt(AABC). □ 

Lemma 7. Let s, t be positive numbers such that s < 1 and t < 1, and let k a positive 
integer. Then 

oo 

JJ(i_ s f) fc > 1-Jfes/(1-*) 

i=0 
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and 



oo 

Y[(l + sf) k < exp(ks/(l-t)). 

i=0 



Proof. The first inequality follows because (1— a)(l— 6) > 1—a—b for a, b G [0, 1], Applying 
this repeatedly, k times for each factor in the product, (1 — s) k {l — st) k ■■■{! — st n ) k > 

l-ks-kst kst n >l-ks J2Zo ? = 1- ks/(l - t). 

The second inequality follows by taking logs and using the inequality log(l + x) < x: 

log((l + s) k ■ ■ ■ (1 + st n f) = k log(l + s) + k log(l + st) + ■ ■ ■ + k log(l + st n ) 

< ks + kst-\ h kst n 

oo 

< ksj^t 1 

8=0 

= ks/(l-t). 

□ 

We now explain how to choose S for the main algorithm. We set it to be 

min{minalt(T) : T G %} 

5 = uw ' (8) 

The minimum altitudes in this definition are measured before any collapse-node operations 
begin. This choice of 5 makes all the theorems work but leads to poorer aspect ratio (by 
a constant factor) than seems necessary. So instead, PINW 1.0 chooses 5 dynamically 
based on the singular values of the affine transformations that could be applied during 
collapse-node operations. This heuristic seems to work well in practice. 

The following theorem bounds the effect of all collapse-node operations, both direct 
and indirect. 

Theorem 5. Let T be a tile in the hierarchy generated by PINW, and let A be the 
composition of all the affine tranformations applied directly to vertices of T and indirectly 
to those vertices via ancestors in the hierarchy. Let a = minalt(T) (prior to any node 
movement). Let I be a line segment lying in T. Assume 5 is chosen according to (0). 
Then length (A(l))/ length(i) lies between 

I 5 V A_ 0.97255\ 3 ^ 0.9725 2 5\ 3 



0.75ay V 0.75a J \ 0.75a 
and 

( 5 \ 3 / 0.9725A 3 A 0.9725 2 <5 

V + 0.75a ) V + 0.75a ) 



0.75a 
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Proof. The proof of this theorem is by induction. The induction base is that for a root 
tile, there are no collapse-node operations so A(l) = I. For a nonroot tile T, let T' be its 
parent triangle. 

By the induction hypothesis, the total distortion of a segment in T' prior to the 
processing of the vertices created within T' is between 

5 V ( 0.9725 A 3 / _ 0.9725 2 5 x 



0.75a'/ V 0.75a' / V 0.75a 
and 

5 \ 3 / 0.97255\ 3 / 0.9725 2 A 3 



l + ■ l + -r^~ • 1 + 



0.75a'/ V' ' 0- 75a ' / V °- 75a 
where a' = minalt(T') (with minalt measured prior to any node movement). Referring to 
Figure 0] and regarding T' = ABC and T as one of I, II, III, IV or V, we consider next 
the direct displacements of D, E, F due to collapse-node operations. 

By Lemma EJ 0.9725a' > a. Thus, for tile T, the distortion prior to the three direct 
displacements of D, E, F is bounded between 

0.97255\ 3 ^ 0.9725 2 ^ y ^ 0.9725 3 ^ y ^ 



0.75a J V 0.75a J V 0.75a 



and 

/ . 0.972V ^ 



0.9725A 3 / 0.9725 2 5V / 0.9725 3 5 



By Lemma [3 with k = 3, s = 6/(0. 75a) and t = 0.9725, infinite product is greater 
than or equal to 1 — 35/(0. 75a (1 — 0.9725)) which simplifies to 1 — 1465 /a. If we assume 
that 5 satisfies (jHJ), then this quantity is greater than 0.9. We now apply the three 
collapse operations of D, E, F to T, say in this order. (Not all three necessarily affect 
T; for example, if T is I in the figure, then moving F does not affect T.) To compute 
the distortion of I requires knowledge of the minimum altitude of T at the point of the 
algorithm when the collapse-node operation is applied. However, because the distortion so 
far is greater than 0.9, we know that the altitude at this step is at least 0.9a for movement 
of D, which is of size 5. Therefore, the movement of D applies a new distortion between 
1 — 5/(0. 9a) and 1 + 5/(0. 9a) by Lemma |U By (jHJ), this quantity is bounded between 
0.95 and 1.05. Therefore, the minimum altitude of T when the collapse-node operation 
for E is applied is at least 0.9a ■ 0.95 = 0.855a. Thus, the collapse-node operation on 
E applies another distortion between 1 — 5/(0. 855a) and 1 + 5/(0. 855a). Again, by (jHJ), 
this quantity is bounded between 0.95 and 1.05. So after the collapse-node operation 
on F, the minimum altitude of T is at least 0.855a ■ 0.95 = 0.812a. Combining these 
three distortions with the distortions from higher-level collapse-node operations given by 
© and (jlUJl shows that the minimum and maximum distortion of a segment after the 
collapse-node operations involving T and its ancestors lies between 

(1 - 5/(0.9a))(l - 5/(0.855a))(l - 5/(0.812a)) 

;1 _ 0.97255 \ 3 / _ 0.9725 2 5 \ 3 / _ 0.9725 3 5 x 3 
0.75a J 'I 0.75a / ' V 0.75a 



25 



and 



(1 + 5/(0.9a))(l + 5/(0.855a))(l + 5/(.812a)) 

/ 0.9725A 3 / 0.9725 2 A 3 / 0.9725 3 A 3 
V + 0.75a J V + 0.75a J \ + 0.75a J 

We can underestimate the first factor and overestimate the second by replacing 0.9, 0.855 
and 0.812 all with 0.75. This proves the theorem. □ 

We now consider Properties 1 and 2 in Section El Since the minimum altitude of a 
triangle is less than or equal to its shortest side length, and since the minimum altitude 
of any tile decreases by at most 0.75, the previous result shows that 5 is sufficiently small 
so that no two nodes can be collapsed to the same node, and no node can have more than 
one choice of where it should be collapsed. 

Furthermore, when we are finished with collapse-node operations, all hanging nodes 
are at least 5 apart and at least 5 from corners. Again, this is because the shortest side 
length is bounded below by the smallest altitude, and the smallest altitude is bounded 
below by a large constant multiple of 5. 

We now consider the aspect ratio of the triangles in the mesh produced by PINW. 
We define the aspect ratio of a triangle to be the square of the longest side length of the 
triangle divided by its area. Since the area is half the product of the longest side length 
and the minimum altitude, an equivalent definition is twice the longest side length over 
the minimum altitude. 

The following lemma gives another characterization of aspect ratio equivalent up to a 
constant factor as well as a useful property of aspect ratios. 

Lemma 8. Let T be a triangle with aspect ratio a. 

(a) Let 9 be the minimum angle ofT. Then there exists two universal constants C\,C2 
such that c\a <l/6< C2a. 

(b) Let l\,l2 be two distinct side lengths ofT. Then lyjli < a/2. 

Proof. Because this lemma is well-known (see, e.g., jl]), we omit the full proof. For (a), 
let H be the length of the longest edge of T. The proof of (a) follows by noting that 
there exists a right triangle one of whose legs has length H and one of whose angles is 
9 that contains T. On the other hand, the same right triangle contracted by a factor of 
1/2 is contained in T. For (b), we observe that I1I2 > 2area(T), i.e., I2 > 2area(T)//i = 
2/i(area(T)//2) > 2l x /a. □ 

The first step of PINW, which performs a preliminary triangulation of Q using Tri- 
angle, outputs triangles that have their aspect ratios bounded above. The reason is that 
Triangle is a guaranteed-quality mesh generation algorithm that will put sharp angles into 
its output only when the input polygon has very sharp angles. Thus, the small angles 
of all the initial triangles have a lower bound. (The reciprocal of the smallest angle of 
a triangle is within a constant factor of the aspect ratio definition given in the previous 
paragraph.) The operation of subdividing at in-centers done to obtain T from Trian- 
gle's output does not increase the longest side length, and reduces the area by at most a 
constant factor. Hence the triangles in % still have bounded aspect ratio. 
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Next, we consider the tiles in 7}, that is, the leaf tiles. Each of these is similar to a 
root tile or its conjugate. In a preliminary step, we ensured that c — a is bounded below 
for all conjugates of root tiles. Therefore, the leaf tiles all have bounded aspect ratio. 

In more detail, the smallest angle of each conjugate leaf tile is either a, where a is the 
smallest angle of a root tile, or is c — a, where a is the smallest and c is the largest angle of 
a root tile. But we have ensured that c — a > .4 by our preliminary splitting rule. Thus, 
if the smallest angle of a conjugate tile is c — a, this means that the conjugate tile has a 
universal upper bound on its aspect ratio. 

Now, we consider the effect of collapse-node operations. 

Lemma 9. After all collapse-node operations are complete, the aspect ratio of any leaf 
tile has increased (compared to its value prior to all collapse-node operations) by at most 
a factor of 1.22. 

Proof. As explained in the proof of Theorem the smallest distortion due to all collapse- 
node operations for any segment in any leaf tile is 0.90 or greater. Pick a tile T and 
let a be the initial altitude. Applying the second part of Lemma [3 to the bound in the 
theorem shows that the maximum distortion of T is exp(3<5/(0.9a(l — 0.9725))), which by 
(JBJ) is at most 1.09. Since the aspect ratio is the twice the longest side length divided by 
the minimum altitude, and the longest side went up by at most 1.09 while the minimum 
altitude changed by a factor at least 0.90, the new aspect ratio is bounded by 1.22 times 
the old. □ 

For this theorem and the remainder of the section, let R\ denote the maximum aspect 
ratio among root tiles and their conjugates. As noted above, because of the properties 
of Triangle, R\ is bounded above by a constant multiple of the reciprocal of the sharpest 
angle of Q. 

Theorem 6. Assume that no root tile is in Tf (i.e., each triangle in % is split at least 
once by the PINW subdivision procedure) . Then, prior to collapse-node operations, the 
maximum value of the minimum altitude among all leaf tiles is no more than cR\ times 
the minimum value of the minimum altitude among all leaf tiles, where c is a universal 
constant. 

Proof. Recall that the tile selected for splitting at any given step is the one with the 
maximum minimum altitude. Thus, when the subdivision procedure terminates, the tile 
at the top of the heap will be the leaf tile with the maximum minimum altitude among 
all leaf tiles. Say this tile is T and its minimum altitude is a. Now consider any other 
leaf tile T' e Tf. Because of the assumption that no tile from T is a leaf tile, this tile 
T' must have arisen from a subdivision of some other tile T". Because of the heap order, 
the minimum altitude a" of T" exceeds a. By Lemma |Hl this means that the minimum 
altitude a 1 of T' is at least min(0.0044, sin a)o/'. Note that sin a > c/Ri since a is an angle 
of a root tile. Thus, a' > (c/R^a" > {c/R x )a. □ 

We now come to the main result for this section about the aspect ratio of the triangles 
generated by PINW. 
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Theorem 7. Each triangle in the simplicial mesh output by PINW has aspect ratio at 
most cR\, where c is a universal constant and R\ was defined above to be the largest aspect 
ratio among root tiles. 

Proof. Let T be a leaf tile. Let a be its minimum altitude and M its longest side length 
prior to any collapse node operations. As already observed in the proof of Theorem at 
the end of collapse-node operations, its minimum altitude is at least 0.9a and its longest 
side length at most 1.09M. The edges of T contain hanging nodes. The distance between 
any pair of adjacent hanging nodes or between a hanging node and corner is at least 5. 
This is because the shortest side length of any leaf tile is a sizable constant multiple of 5, 
so no leaf tile edge could ever shrink below 5. 

Recall from (JBJ) that 5 is a constant multiple of a m i n , the minimum altitude among 
all leaf tiles. By Theorem EJ this implies 5 > ca/Ri, where c is a universal constant and 
a = minalt(T). 

Now, let r be a triangle output by the Delaunay triangulation of T and consider the 
sharpest angle of r. Let e = Viv 2 be edge of r opposite the sharpest angle. There are two 
cases: either vi,v 2 lie on the same edge of T (i.e., they are consecutive hanging nodes or 
a hanging node and a corner node), or they are on different edges. 

Start with first case. Let e = viv 2 be the edge of r lying on an edge of T. For the rest 
of this case, let T = A ABC such that V\V 2 C BC and such that the order of these vertices 
is B, V\,V2, C . As mentioned above, |e| > cajR\. Let p be the vertex of r opposite e. By 
definition of the Delaunay triangulation, p is the first vertex hit by an expanding circle that 
contains both endpoints of e. This circle, if expanded further, will eventually hit A, the 
vertex of T opposite the edge of T containing e. Either ZAv\v 2 is acute or ZAv 2 v\ is acute; 
without loss of generality, assume the former. The angle Zv\Av 2 is at least c\e\/(RiM) by 
the law of sines: sin Zv\Av 2 = sin ZAviV 2 \e\/\Av 2 \ > sin ZAviV 2 \e\/(1.09M). Meanwhile, 
sin ZAviv 2 is bounded below by c/Ri since ZAv\V 2 is bounded below by ZABC but is less 
than 7r/2. Thus, sinZv^v-z > c\e\/(RxM). Next, M < aRi so sin Zv x Av 2 > c\e\/(aR\). 
Finally, \e\/a > cjR\ as noted in the previous paragraph. Therefore, sin Zv\Av 2 > c/Rf. 
This is the angle formed by viAv 2 . The actual angle of r opposite e is vipv 2 . But this 
angle is greater than or equal to v\Av 2 , since the expanding Delaunay circle encounters p 
before it encounters A (or else p = A). 

Next, let us consider the case that v\ and v 2 , the endpoints of the edge of r opposite 
its sharpest angle, do not lie on the same edge of T. Let the three vertices of T be A, B, C 
and let A be the T-vertex that is the common endpoint of the two T edges that contain 
v i and v 2 respectively, while we let B be the T-vertex such that AB contains v i and we 
let C be the T-vertex such that AC contains v 2 . Consider Zv\Bv 2 . Since A, v\ and B are 
collinear, this angle is equal to ZABv 2 . By the law of sines applied to AABv 2 , we have 
sin ZABv 2 = \Av 2 \sin ZA/\Bv 2 \. Now we apply the following inequalities: \Av 2 \ > 5, 
sinZA > c/Ri, and \Bv 2 \ < M to conclude that sinZABv 2 = sinZv\Bv 2 > c5/(RiM). 
This was the same inequality derived in the previous case, and yields the conclusion that 
sin Zv\Bv 2 > c/R\. Arguing again as in the previous case, the actual Delaunay triangle 
containing v\v 2 may not have B as its third vertex, but if it has any other vertex w, then 
Zv\wv 2 is greater than Zv\Bv 2 by considering the expanding circle property. □ 
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8 Isoperimetry of the final mesh 



We have already proved in Theorem 01 that the tiling of a triangle by our generalized 
pinwheel subdivision has the isoperimetric property. It is straightforward to extend this 
result to the collection of all leaf tiles. 

Theorem 8. Let %, 7[, . . . , be the sequence of tilings of Q generated by the PINW algo- 
rithm as follows. For each n, T n is the set of leaf tiles of Q generated by PINW when the 
user- specified size requirement is 5 n such that 5 n — > as n — > oo. Then for any distinct 
points P, Q such that P,Q G Skel(Tfc) for some k, we have 

lim dist (P,Q) = \\P-Q\\n. 

n^oo Skel(T n ) 
n>k 



Proof. This theorem follows from Theorem 0] and uses the same proof technique. Let II be 
the geodesic path from P to Q of length \\P — Q\\q. Since Q is a polygon, II is composed 
of a finite number of line segments. For each tile Tj in 7& that meets II, consider the small 
segment PiQi that is Tj n II. Then we use Theorem El to argue that this small segment 
PiQi can be approximated arbitrarily accurately. □ 

This theorem can now be extended to the final mesh output by PINW by analyzing 
the effect of collapse-node operations on the isoperimetric number. (The Delaunay op- 
erations do not disturb the isoperimetry result, since adding edges could only make the 
isoperimetric number decrease.) 

The definition of isoperimetry implicit in Theorems|3]and|B]is not suitable for analyzing 
the output of PINW because the meshes produced by PINW are not refinements of their 
predecessors as the mesh size decreases. This is because the collapse-node operations 
move nodes differently depending on the size of the leaf tiles. 

Therefore, we use the following definition. An infinite sequence of simplicial meshes 
M.2) ■ ■ ■ for a domain Q has the isoperimetric property if for each A4i there is a subset 
Li of its vertices such that the following two properties hold. First, Tj is asymptotically 
dense in Q as i — > oo, i.e., for any e > 0, there is an / such that for any x G Q and any 
i > I, there exists a v G Li such that \\x — v\\ < e. Second, 

j dist SkcKMi) (x,y) \ 

lim sup < r r : x, y G 1^; x ^ y > = 1. 

»-> 00 I \\x-y\\n J 

Theorem 9. The family of meshes produced by the PINW algorithm has the isoperimetry 
property described in the previous paragraph. 

Proof. To show that PINW has this property, take a sequence of q's tending to zero. 
For each i, let % be a generalized pinwheel subdivision of Q such that each leaf cell has 
diameter less than e,/2. Then let T? be a further subdivision of % such that for any two 
distinct vertices of %, dist Skel (-7-/)(£ — y) < (1 + ej/4)||x — y\\Q. The existence of such an 
T? is established by Theorem |H1 Let a' be the minimum altitude among leaf tiles in 7^'. 
Next, further refine T[ to yield a tiling T" with the property that when 5 is defined by 
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(JBJ for T", (i.e., the % appearing in (JHJ) pertains to %"), then this S is sufficiently small 
so that 

exp(35/(0.75a'(l - 0.9725))) < 1 + e,/4 (11) 

and 

1 - ei/4 < 1 - 35/(0.75a'(l - 0.9725)). (12) 

Now finally, take J\A% to be the simplicial mesh output by PINW based on T", and take 
Li to be the set of nodes of M.i that are displaced copies of the nodes of %. 

First, we have to show that Lj defined in this manner is asymptotically dense. The 
nodes of Lj are the same as the nodes of % after small displacements. Since every cell of 
% has diameter less than ej/2, this means that any point x G Q is distance at most 
from a vertex of %. The vertices of Lj are slightly displaced, but no distance d decreases 
below 0.9d nor increases to more than 1.09c?. Therefore, for any igO the perturbed set 
Li contains a point v within distance 1.09 • ej/2 < q of x. 

Let x, y be two distinct points in Lj. The next task is to show that distgk e i(A^ i )(^, y) < 
1 1 ar — y||n(l + e i)- Let x , yo be the positions of x, y in 7?' prior to all distortions caused by 
collapse-node operations. Note that x , y are vertices of 7? and also of % by construction. 
Therefore, by construction of T/, there is a path P in Skel(7~') connecting x and y such 
that length(P ) < ||x — yo||n • (1 + q/4). Let the segments of P be ei, e 2 , . . . , e r . Let 
the image of Pq after all the collapse-node operations with their attendant distortions are 
applied be P, and the images of e\, . . . , e r be fx, ... , f r . Recall that the distortions that 
affect a node v of a tile T in the hierarchy are those distortions associated with T and 
its ancestor tiles, but descendant tiles cannot move T. Therefore, by Theorem all of 
the quantities \\x — y\\/\\xo — yo\\u, length (/j)/ length(ej) , and length(P)/length(P ) he 
between 

5 \ 3 _ 0.9725A 3 / 0.9725 2 (T 



0.75a7 V 0.75a' J \ 0.75a 
and 

5 \ 3 ( 0.97255\ 3 / 0.9725 2 5 
1 + 777^ 1 • I 1 + . J ■ I 1 + 



0.75a'y V °- 75 «' / V °- 75 «' 
where the 5 in this formula is given by (jHJ) associated with T" . By Lemma [7[ this interval 
is bracketed by 

1 



0.75a'(l - 0.9725) 
and 

exp(3<5/(0.75a / (l - 0.9725))). 
Then by (fTTj) and (|T2|) . this interval is bracketed by [1 — ej/4, 1 + Cj/4]. Thus, 

dist (x,y) < (l + ei/4) dist (ac ,2/o) 

Skel(Mi) Skel(T/) 

< (1 + ei /4) 2 ||x -y \\n 

< (l + e i /4) 2 ||a;-y||n/(l-e l /4). 

Note that (1 + ei/4) 2 / (1 — e*/4) < 1 + as long as 6j < 1/2. Thus, we have shown that 
for all x,y E L h 6.ist^ Mi) {x, y) < (1 + e<)||a: - y|| n - □ 
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9 Conclusions 



We believe that this generalization of pinwheel tiling to meshing polygonal regions would 
aid in modeling arbitrary crack paths more accurately than the current meshing tech- 
niques. Also, this kind of substitutive mechanism for subdivision makes it easy for adap- 
tive meshing. This work raises a number of interesting directions for future research. 
Among them are the following: 

1. The transformation of the tiling to the mesh had the effect of increasing the aspect 
ratio significantly. Is there a better way to carry out this transformation to reduce 
the impact on aspect ratio? 

2. The convergence rate of the isoperimetric number of the pinwheel tiling to 1, which 
was not analyzed here, is known to be extremely slow even in the case of the 1:2 
tiling. Is there another approach to isoperimetry that converges faster? 

3. Consider a mesh generated by placing random points in the domain under con- 
sideration and joining them with a Delaunay triangulation. Is there a limiting 
isoperimetric number for this family of meshes (with high probability)? 

4. Another way to construct a mesh of an arbitrary polygon with limiting isoperimetric 
number equal to 1 is to use the 1:2 pinwheel subdivision for every coarse triangle 
after subjecting it to a (potentially large) afline transformation. This approach is 
simpler in certain respects than PINW. For example, the collapse-node operations 
for this algorithm need to be done only at the boundaries of the coarse triangles. The 
difficulty with this approach is that it spoils the "statistical rotational invariance" of 
the pinwheel tiling. The statistical rotational invariance property states that the set 
of possible directions is covered at a uniform rate as subdivision proceeds. We are 
unclear whether statistical rotational invariance is important for cohesive interface 
modeling. We suspect that our construction of generalized pinwheels has statistical 
rotational invariance but have no proof of this. 

5. Can any of this work be extended to three dimensions? 
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